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Abstract 

We study the evolution under their self-gravity of particles evolving from infinite "shuffled lattices" 
initial conditions. We focus here specifically on the comparison between the evolution of such a 
system and that of "daughter" coarse-grained particle distributions. These are sparser (i.e. lower 
density) particle distributions, defined by a simple coarse-graining procedure, which share the same 
large-scale mass fluctuations. We consider both the case that such coarse-grainings are performed 
(i) on the initial conditions, and (ii) at a finite time with a specific additional prescription. In 
numerical simulations we observe that, to a first approximation, these coarse-grainings represent 
well the evolution of the two-point correlation properties over a significant range of scales. We 
note, in particular, that the form of the two-point correlation function in the original system, 
when it is evolving in the asymptotic "self-similar" regime, may be reproduced well in a daughter 
coarse-grained system in which the dynamics are still dominated by two-body (nearest neighbor) 
interactions. This provides a simple physical description of the origin of the form of part of the 
asymptotic non-linear correlation function. Using analytical results on the early time evolution of 
these systems, however, we show that small observed differences between the evolved system and its 
coarse-grainings at the initial time will in fact diverge as the ratio of the coarse-graining scale to the 
original inter-particle distance increases. The second coarse-graining studied, performed at a finite 
time in a specified manner, circumvents this problem. It also makes more physically transparent 
why gravitational dynamics from these initial conditions tends toward a "self-similar" evolution. 
We finally discuss the precise definition of a limit in which a continuum (specifically Vlasov-like) 
description of the observed linear and non-linear evolution should be applicable. This requires the 
introduction of an additional intrinsic length scale (e.g. a physical smoothing in the force at small 
scales), which is kept fixed as the particle density diverges. In this limit the different coarse-grainings 
are equivalent and leave the evolution of the "mother" system invariant. 

PACS numbers: Pacs: 05.40.-a, 95.30.Sf 



I. INTRODUCTION 

The dynamics of infinite self-gravitating systems of 
point particles starting from quasi-uniform initial con- 
ditions is a straightforward, well posed, but essentially 
unsolved problem of out of equilibrium statistical me- 
chanics. While numerical simulation of such systems has 
developed impressively in scale and sophistication in the 
last few decades — mostly driven by the relevance of the 
problem in cosmology — analytic understanding of the 
evolution of the clustering observed remains very limited. 
The work reported in this article is a continuation of a 
study of this problem reported by us recently M. We 



study a simplified version of the problem posed in cos- 
mology, considering pure Newtonian gravity in a static 
cuclidcan universe and a very simple class of "shuffled 
lattice" initial conditions (IC). 

In Q we have presented in detail the results of a numer- 
ical study of this system. We have found that, as soon as 
strong positive two-point correlations begin to develop, 
they are characterized by a simple spatio-temporal scal- 
ing relation, in which the strong spatial correlations at 
any time may be inferred from those at a precedent time 
by a simple rescaling of the spatial variables. Further 
the function of time specifying this dynamical scaling of 
the characteristic spatial scale tends, after a transient 
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period, to a form which can be determined solely from 
the linearized fluid equations describing the evolution of 
the small fluctuations at large scales. This asymptotic 
dynamical scaling behavior is completely analogous to 
that which has been observed in cosmological simula- 
tions, starting from more complicated IC (correlated per- 
turbations of lattices) and in an expanding universe, and 
referred to in this latter context as "self-similarity" . The 
amplitude of the two-point correlation function is also 
observed in this asymptotic regime to be essentially inde- 
pendent of the mean particle density over a large range of 
scale. This suggests that the full non-linear dynamics of 
these simulations, once this asymptotic scaling behavior 
is approached, might be well approximated solely within 
the framework of a continuous fluid-like description of 
the system. We have found, however, that the form of 
the two-point correlation function (i.e. its dependence on 
separation rather than its temporal scaling behavior) at- 
tained in this asymptotic scaling regime is already very 
well approximated, at least for modest amplitudes, by 
that at early times where the dynamics at these scales 
is explicitly particle-like in nature. Indeed the correla- 
tion function which develops at these times is very well 
approximated by assuming that the dynamics is that of 
particles interacting only with their few nearest neighbor 
particles. In other words, while the spatio-temporal scal- 
ing law seems to be determined by only the large scale 
fluid linear dynamics, the form of the scaling function 
looks to depend on the discrete aggregation dynamics 
between elementary objects. This suggests that the dis- 
creteness of the system may be an essential ingredient 
even to understand its long-time behavior. 

In this paper we study further the dynamical evolu- 
tion in this same class of models (i.e. with "shuffled 
lattice" IC in a static universe). Specifically we focus 
here on the comparison between the evolution of a given 
such IC and the evolution of other point distributions de- 
fined from them by a simple coarse-graining procedure: 
at a chosen time the mass in coarse-graining cells, de- 
fined by a lattice, is aggregated into a single mass; this 
new point mass is attributed the position and velocity of 
the center of mass of the mass in the cell. This is a proce- 
dure which evidently modifies the system up to a certain 
scale — that characteristic of the coarse-graining — but 
leaves the mass and velocity fluctuations at larger scales 
essentially unperturbed. The comparison of the origi- 
nal and coarse-grained (CG) system allows one to infer 
useful physical information about the dynamics. Specif- 
ically it allows one to study the role of fluctuations at 
different scales in the system, and the degree to which 
a discrete or continuous description of the dynamics can 
be appropriate. We will see that the study will allow us 
to understand better the physics underlying the asymp- 
totic "self-similar" evolution, and also to address more 
precisely the question of the continuum limit of the sys- 
tem. These latter considerations are of interest, as we 
will explain in the conclusions, in the context of simu- 
lations in cosmology. Such simulations can be seen as 



particle coarse-grainings in the sense defined here of the 
"true" self-gravitating matter described by cosmological 
models 1 . A clear understanding of the continuum limit is 
important as it should correspond to the theoretical cos- 
mological model. Discreteness effects in these simulations 
are the differences between the CG and this underlying 
continuum model. 

In the introduction to |jj we have given a brief review of 
the relevant literature, on both self-gravitating systems 
and other long-range interacting systems. The notion 
of "particle coarse-graining" , which is the central one in 
this paper, is simply a specific kind of coarse graining 
procedure introduced in statistical physics (see, e.g.|||). 
It is evidently not unnatural to consider such schemes in 
approaching the dynamics of systems manifesting scal- 
ing behavior, as has been done widely in many different 
systems in condensed matter physics with well known re- 
sults. In the context of self-gravitating systems various 
authors have studied self-gravitating dynamics applying 
such concepts (see, e.g., ||, |J [|, H, @)- These works, how- 
ever, study theoretical descriptions of these systems in a 
continuum limit. In this work we consider, instead, the 
fully discrete systems, numerically simulated, applying 
coarse-graining procedures which relate one such discrete 
system to another. Indeed one of the questions we ad- 
dress in this work is, as mentioned above, the validity of a 
continuum description of the dynamics. One of our find- 
ings is that the coarse-grained system, in a phase in which 
its dynamics are manifestly discrete, already traces well 
the behavior of the original system. The only previous 
work we are aware of which has actually applied a discrete 
"renormalization" procedure is that of ||, |9|, which at- 
tempts to determine the mass distribution which results 
from power law initial conditions using such a scheme 
under the assumption that the evolution is self-similar. 
Numerical studies of discreteness effects in cosmological 
N body simulations (e.g. [^0)) also implicitly consider 
the effect of different discrete coarse-grainings, albeit not 
defined as such in the simple way considered here. In- 
deed, as mentioned above, one of the goals of this work 
is to clarify problems in this context by a well controlled 
study of much simpler systems. 

The paper is organized as follows. In the next section 
we summarize briefly the essential definitions and the re- 
sults of the preceding paper [[[) . We also summarize rele- 
vant results derived recently by the present authors and 
B. Marcos fH |l| on a perturbative description of the 
early times dynamics of self-gravitating points perturbed 
off a lattice. We underline the slight differences in the 
systems studied here compared to that in jij — a differ- 
ent probability distribution function for the initial "shuf- 
fling" and also a different ascription of initial velocities. 



For typical "cold dark matter" models there is a factor of the or- 
der of 10 70 between the microscopic number density of the phys- 
ical particles and the number density of the "macro-particles" 
used in the largest simulations performed up to now. 
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To explain this latter point we include a short discus- 
sion of the "Zeldovich approximation" to the evolution 
of a pressure-less self-gravitating fluid (This approxima- 
tion is also important for the discussion in Sec. |V| of the 



paper). In Sec. Ill we define precisely and discuss the 
construction of the coarse-grained particle distributions 
from a given initial particle system. In the following sec- 
tion we report our study of the comparison between a 
given SL and a set of CG of the initial conditions, while 
m Sec.0 we do the same for a set of CG defined at a given 
scale at finite time in the course of the evolution. In the 
final section we summarize our results and conclusions, 
and discuss what we have learned about the dynamics 
of these systems. Specifically we discuss the question of 
the appropriateness of a continuum description of the ob- 
served non-linear dynamics, and we indicate some paths 
for future research on this question. 



II. BASIC DEFINITIONS AND BACKGROUND 
RESULTS 

In this section we give the essential definitions and 
background for the paper. Firstly we summarize the ba- 
sic definitions of the initial conditions studied and of the 
quantities used to characterize the systems. Except for a 
few minor points this discussion essentially summarizes 
an analogous one in where the interested reader can 
find further details. One minor difference is that wc con- 
sider here the analysis also of particle distributions in 
which the masses of all points are not equal. We then 
summarize the essential results of jlj , and then also those 
of which we also make use of in this paper. Finally, 
we discuss briefly the so-called "Zeldovich approxima- 
tion" which is an important perturbative approximation 
to the evolution of self-gravitating systems in the fluid 
limit, and therefore an approximation to the evolution 
of the particle system in a certain regime (long wave- 
lengths and small amplitude perturbations). We will use 
it specifically, as it is commonly done in cosmological sim- 
ulations, in the choice of initial conditions for the particle 
velocities, and also at various points in our analysis. 



A. Definition of a Shuffled Lattice 

We call shuffled lattice (SL) the particle distribution 
obtained by applying a random displacement indepen- 
dently to each particle on a simple cubic lattice. Thus 
the position of a particle, initially on the lattice site R, is 
given by x(R) — R + u(R) where the vectors u(R) are 
random and specified by a probability density functional 



V[{u(R)}]=l[ P [u(R) 



(1) 



sity function (PDF) for the displacement of a single par- 
ticle. 

In the numerical simulations reported in jjj we took 
a simple "box" form for the PDF. Here instead, for rea- 
sons which we will explain below, we will take p(u) to be 
simply a mono-variate Gaussian: 



p(u) = p{u) = 



1 



27TCT 2 



3/2 



o.xp(-^). (2) 



where u is the modulus of u. As discussed in Q we 
expect the choice of the precise form of the PDF (for 
any simple functional form) to have little effect on the 
dynamics beyond some early time transients. 

We use the term shuffling parameter, denoted A, for 
the rms deviation of it, i.e., 



A^ 



u 2 p(u) d 3 u = 3cr 2 



(3) 



We will usually use this quantity expressed in terms of 
the lattice spacing, which we denote by £. We thus de- 
fine the normalized shuffling parameter S = A/£. This 
is convenient because it is actually the value of S alone 
which characterizes the gravitational dynamics of an in- 
finite SL M: any two SL with the same S but different 
t (and A) are equivalent up to a scale transformation 
which is irrelevant for (scale-free) gravity. We note that 
in the limit <5 — s> the particle distribution remains a per- 
fect cubic lattice, while for S — > oo it becomes a Poisson 
distribution. 



B. Characterization of particle distributions 

We now discuss briefly the different statistical quanti- 
ties which we will use to characterize the SL initial condi- 
tions and the evolved particle distributions. For further 
details, we refer the reader to [§ (D§ @. 

For any distribution of N particles in a volume V , we 
can define the microscopic mass density as 



p(x) 



N 



i=l 



mi5v(x 



(4) 



where Xi S V is the position of the ith particle, of mass 
mi, and i5d is the Dirac delta function. For the case of 
infinite systems with a well defined non-zero mean den- 
sity po, which is that we consider here, it is convenient 
to define the density contrast S p (x) — [p(x) — po}/ Po- 

The two-point correlation properties can then be char- 
acterized by the reduced two-point (density-density) cor- 
relation function 



R 



Ur) = (S P {r + x)5 p {x)) 



(5) 



The statistical properties of the particle distribution are 
thus completely determined by p(u) , the probability den- 



where (. . .) is an ensemble average. In our simulations we 
treat N particles in a box of side L = N 1 / Z l with periodic 
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boundary conditions. We can therefore write the density 
contrast as a Fourier series: 



Sp(x) = -3 ^2 exp(ife • x)Sp(k) 



(6) 



with k € {(2?r / L)n\n S Z 3 }, and 



<5 p (fc) = / exp(— ik ■ x)S p (x) 
II 3 



1 

Po 



exp(— ife • Xj) — 8^{Lk/2-K, 0) : 



(7) 



where 5k is the Kronecker delta function. The power 
spectrum 2 (PS) may then be defined (see, e.g. |14|) as 



P(k) 



1 

u 



\S P (k)\ : 



(8) 



For statistically homogeneous (i.e. statistically transla- 
tional invariant) particle distributions 3 , the PS and £(r) 
are a Fourier conjugate pair and therefore contains the 
same information. 

We estimate the PS as a function of the modulus k = 
I A; I, using simply the expression given in Eq. (Q), averaged 
over all k in a shell in k space. We therefore denote it 
from now on as P{k). For the real space correlations we 
will consider estimating £ p (r) as defined above only in 
the cases where the points of the distribution considered 
have equal mass. Just as for the PS we will evaluate it by 
performing an average over all r in a radial shell. In this 
case we will denote it as £,(r), and it can be written |jl4), 
for r 7^ 0, as 



n 



(9) 



where (n{r)) is the radial conditional mean number den- 
sity (i.e. the mean number density of particles at a dis- 
tance r from an occupied point), and no is the mean (un- 
conditioned) number density. This makes its direct esti- 
mation in real space very simple, by simple pair-counting 
algorithms (see W for the explicit description of the al- 
gorithm we choose). 

Besides these two quantities we will consider also, as in 
jjj, one other one: the nearest neighbor (NN) distribution 
w(r), which is simply the probability density function for 



2 As in [y we use the terminology and definitions current in cos- 
mology. In solid state physics the PS is usually called the struc- 
ture factor, and normalized typically so that it is dimcnsionlcss 
rather than having units of volume. 

3 In a lattice, the ensemble average is defined over the set of lattices 
translated by an arbitrary vector in the lattice cell. In the case 
of the SL we have both this average and, in addition, that over 
all realizations of V[{u(R)}]. 



a particle to have its NN at a distance r. As in [Q it is 
estimated in the evident way by direct counting of NN 
separations. The usefulness of this quantity here is that 
it allows one to infer |]l5| very useful information about 
the nature of the correlations at small scales, and thus 
indirectly about the dynamics responsible for them: if 
we neglect all but the two-point correlations one has the 
relation (TJ] 

u(r) = (l - uj{s) ds^j ■ 47rr 2 n [l + £(r)] . (10) 

We will use this relation (as in jl], |l5| ) to probe the degree 
to which interactions with NN particles are responsible 
at early times in simulations for the development of the 
observed correlations. 



C. Gravitational clustering in a Shuffled Lattice 

It is straightforward [jfj], |l4|, to calculate exactly 
the two-point correlation properties of the SL. The PS is 
given by 



P(k) 



l-|p(fc)| ; 
n 



9tt 

L 3 ]T|i5(fc)| 2 Mfc,n^) (11) 



where p(k) is the Fourier transform of the PDF for the 
displacements p(u) (i.e. its characteristic function). For 
the case of a Gaussian PDF as in Eq. (S), with Eq. (^), 
it follows that, for k ^=n=f-, 



P(k) = P(k) = — 
n 



1 — exp 



A 2 fc 2 



At leading order in small k we therefore have 



P(k) = P{k) 



2 1,2 



A 2 fc 
3n n 



(12) 



(13) 



Note that, at large k, one has P(k) — 1/no, which is 
simply the particle shot noise necessarily present in any 
distribution of point-particles. 

In H we have studied the dynamical evolution of self- 
gravitating particles starting from SL initial conditions, 
and zero initial particle velocities. The PDF p(u) used 
for the initial displacements off the lattice is a simple 
"box" form, i.e., with uniform probability inside a cube 
centered on the origin and oriented parallel to the lattice. 
The leading term in the PS Eq. (|l3|) is in fact independent 
of the form 4 of p(n), and thus in this approximation the 
initial conditions are the same. 

The principal results of M are the following: 



4 The only assumption about the PDF is that is has finite variance. 
The qualitatively different case of infinite variance PDF is treated 



also in 



0,0. 
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• Qualitatively the evolution of clustering is very 
similar to that observed in cosmological simula- 
tions (see, e.g., Jl7| ) , which are performed in an 
expanding spatial background and start from lat- 
tices with correlated displacements representing PS 
typically with P(k) oc k n and —3 < n < —1. 
Clustered structures develop initially at very small 
scales (< t) and then progressively at larger and 
larger scales. 

• The evolution is observed to be independent of the 
box size L until the time when the size of the 
non-linear structures approaches the box size itself. 
Thus the results of numerical simulations arc inter- 
preted as representative of the infinite volume limit 
(taken at constant particle density). 

• The PS at small k is amplified as predicted by the 
linearized fluid limit for the system (see, e.g. [pi). 
The k below which this behavior is observed de- 
creases with time, reflecting the propagation of the 
non-linear clustering in real space to larger scales. 

• The temporal evolution of two-point correlation 
function which coincides in this case (of equal mass 
points) with the density-density correlation func- 
tion, is well described from early times by a spatio- 
temporal scaling relation 



fat) 



f {\(t) 



(14) 



i.e., the temporal evolution of £(r, t) is well approx- 
imated by a simple rescaling of the spatial coordi- 
nates 5 . 

After a time At this scaling behavior converges to- 
ward a more specific form which may be written 



£(r,i + At) = £ 



R s (At) 



R s (At) = e , (15) 



where Td yn is the dynamical time of the system, de- 
fined by Tdyn = 1/V4ttG/9o- The constant z, which 
one can consider a sort of dynamical exponent char- 
acterizing the dynamical scaling, is given by 

* = T~~ — ( 16 ) 



where n is the exponent of the power-law PS of the 
initial particle distribution (and therefore z = 2/5 
for the SL). This behavior is completely analogous 
to what is referred to as "self-similarity" in the con- 
text of cosmology (see, e.g., [0). It is the behav- 
ior which follows if one assumes that (i) the char- 
acteristic length scales in the system play no role 



in the dynamics, and (ii) the linearized fluid limit 
correctly describes the evolution of fluctuations at 
low amplitude. It is thus expected to be observed 
specifically starting from IC which are power law 
PS, if the clustering is independent of the lattice 
spacing (which is, in this case, the sole character- 
istic length scale). The results in Q extend the 
range in which such behavior has been reported in 
numerical simulations in cosmology to "bluer" PS 
(and also to the case of a static universe). 

• We have mentioned that the sole relevant param- 
eter for the gravitational dynamics of a SL is the 
normalized shuffling parameter S. In jj| different 
SL with a range of S, but identical large scale am- 
plitude of the PS, are compared. The evolution of 
these IC are all observed to converge to the self- 
similar behavior after a transient which depends 
on S. These results are closely related to those de- 
scribed in Sec. IV below. 



5 In the rest of the paper we will often write § (r) leaving the time 
dependence implicit. 



• In the early phase of evolution, in which non-linear 
clustering develops, but prior to the onset of the 
self-similarity, the observed two-point correlations 
are very well accounted for solely by two body cor- 
relations developing under the influence of NN in- 
teractions. The form of the two-point correlation 
function at this time is also a rather good approxi- 
mation to that in the self-similar regime. Indeed, as 
noted above, the spatio-temporal scaling of the cor- 
relation function is a good approximation well be- 
fore the self-similar regime is attained, i.e., in the 
transient phase we have already a temporal evo- 
lution given approximated well by Eq. Jl4|), but 
with a different temporal dependence of the func- 
tion X(t) to that in the asymptotic regime of dy- 
namical scaling Eq. 

The aim of this paper, as has been discussed in the 
introduction, is to gain more insight into the dynamics 
at work in these simulations by analyzing the effect of 
coarse-graining the system at different times in the course 
of its evolution. 



D. Early time evolution of a perturbed lattice 

In |ll], [l2| we have developed in detail a general formal- 
ism for treating the evolution of self-gravitating particles 
perturbed off an infinite perfect lattice, in the limit that 
these perturbations are such that the relative displace- 
ments of the points are small compared to their separa- 
tions. The treatment is thus valid up to times at which 
this approximation breaks down, and is applicable to any 
set of perturbations, whether uncorrelated (as considered 
here) or correlated (as in cosmological simulations). It 
applies equally well to the case of a static euclidean uni- 
verse or a cosmological expanding one. We will make use 
of these results below and thus summarize them briefly 
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here. We refer to this approximation to the full dynamics 
as particle linear theory (PLT). 

The formalism is essentially the same as that used 
canonically to analyze classical phonons in a crystal in 
solid state physics (see e.g. Jigj] ) . One can write the 
equation of motion for the displacement u(R) of a par- 
ticle, from its lattice site R, as: 



u(R, t) = - E £>(R - R') U (R', t) 



(17) 



R' 



where T> is the dynamical matrix, which is derived by 
linearizing the force on a particle in the displacements 
relative to all other particles. For gravity we have 



2V(R^0) = Gm 



ZV(O) 



R :i 

E 2V( R ) 

R^O 



R 5 



(18) 
(19) 



Defining the discrete Fourier transform on the lattice by 



u(k,t) 



E« 

R 



ik-R 



u(R,<) 



u(R,t) = ^E eik ' Rfi ( k '*) 



(20) 
(21) 



where the sum in k is over the first Brillouin zone, i.e., 
for the simple cubic lattice k £ n(2?r/L), where n is a 

vector of integers each g] — N ^ , N 2 ], the equation of 
motion Eq. (pj) becomes, in reciprocal space, 

u(k, t) = -X>(k)u(k, t) (22) 

where 2?(k), the FT of 2?(R), is a real symmetric 3x3 
matrix for each k. 

The dynamical problem thus reduces to the diagonal- 
isation of 2?(k) for each k in the first Brillouin zone. 
Labeling the three orthonormal eigenvectors e„(k) and 
their eigenvalues w^(k) (n — 1, 2, 3), the evolution of the 
displacement field from t = is given as 6 

u(R, t) = ± E k E„ e„(k) • {u(k, 0) • e„(k) cosh[w n (k)i] 



~[k) 



u(k, 0) • e„(k) sinh[cj„(k)i]}e 4k - R (23) 



Thus given the initial displacements and velocities, the 
dynamical evolution is solved. 

The solution of the diagonalisation problem is numeri- 
cally straightforward. Details of it, and the results for the 
eigenvalues and eigenvectors, are given in 11 . In this pa- 
per the domain of validity of this perturbative treatment 
has also been investigated using numerical simulations. 



Starting from (i) an SL initial condition and (ii) from 
a perfect lattice configuration with correlated perturba- 
tions corresponding to a PS ~ k~ 2 , it is found that the 
treatment traces very well, in both cases, the full evolu- 
tion of the system (i.e. with the full gravitational interac- 
tion) until the average relative displacement approaches 
the lattice spacing. 



E. The Zeldovich Approximation 

Keeping the mean mass density po fixed, and taking 
the limit I — > 0, the eigenvectors and eigenvalues at 
any given fixed k simplify: one obtains a single non-zero 



eigenvalue independent of k, lu 2 (k) = inGpo = r dy ^ , as- 
sociated to a longitudinal mode e(k) = k, and two trans- 
verse modes (i.e. e(k)-k = 0) with zero eigenvalues. If we 
consider the case that the perturbations to the lattice are 
of long wavelength, i.e., that there are perturbations only 
with k£ -C 1, we obtain from Eq. ( |23] ) that there is an 
asymptotic attractor solution for large times (t ^> Td yn ), 
which may be written 



u(R, t) = exp(i/Tdyn) q(R) : 



(24) 



where q(R) is the time independent curl-free (irrota- 
tional) vector field: 

l( R ) = Jf E {( fl ( k > °) ' k ) + ^y„(u(k, 0) • k)} • ke lk R 

1. (25) 
Using Eq. (|22| ) it is easy to show that the gravitational 
field can be written, in the same long-wavelength limit, 
as 



g(R,t)=r dy 2 n u(R,t) 



(26) 



so that it follows that the asymptotic solution can be 
written as 

u(R, t) = exp(t/T dyn )r d 2 yn g(R, 0) , (27) 
u(R, t) = exp(i/T dy „)T dyn g(R, o) , (28) 

i.e., both the displacement and velocity of each particle is 
expressed solely in terms of the gravitational field acting 
on it at the initial time 7 . 

This approximation corresponds to one introduced by 
Zeldovich [ p0[ , for the evolution of fluid elements in a 
pressure-less self-gravitating fluid away from a perfectly 



We give here the formulae for the specific case of a static eu- 
clidean universe. The eeneral result for cosmological back- 
grounds is given also in Jill. 



7 The approximation is thus sometimes described as "ballistic" . 
This, however, is a misnomer as the approximation incorporates 
in fact the collective effect of the motion of all particles on one 
another: the approximation is non-trivial precisely insofar as it 
extends well beyond the regime of a ballistic approximation. 
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uniform state, in an expanding universe. It can be for- 
mally derived as an approximation through a pertur- 
bative treatment of the fluid equations in the Lagrangian 
formalism 8 . 

Note that using the continuity equation for small dis- 
placements u applied to elementary volumes in a con- 
tinuous medium one obtains that the density fluctua- 
tions Sp with respect to the mean density po are given 
by Sp « — poV ■ u. Linear amplification of small den- 
sity fluctuations is thus associated with the growth of 
displacements described by Eq. (^4j), and indeed it was 
on the basis of this observation that Zeldovich proposed 
his ansatz. The power of the Zeldovich Approximation 
(ZA) is that its domain of validity extends well beyond 
that of linear (Eulerian) perturbation theory of the self- 
gravitating fluid. In fact it extends up to "shell-crossing" , 
when fluid elements contract in one direction to produce 
density singularities. In terms of our particle treatment 
this shell-crossing corresponds to the approach of parti- 
cles to one another. 

In cosmological simulations the ZA — in the form of 
Eq. ( p7j ) adapted to an expanding universe — is used to 
set up initial conditions: it gives a prescription for both 
the displacements and velocities of the particles off the 
lattice (considered as the centers of fluid elements) once 
the input PS is given, as any realization of the latter 
gives the gravitational field through the Poisson equa- 
tion. Our initial conditions are not equivalent to this as 
the displacement field applied to the lattice at the ini- 
tial time is not, in general, curl- free as in the ZA (since 
the gravitational field is the divergence of a scalar field). 
We can, however, use the ZA to determine the velocities 
since Eqs. (|27]) and ([28|) imply 

u(R,0) = ^M. (29) 

7"dyn 

This is the prescription we will adopt here (rather than 
the zero initial velocities of Q). It is a very natural 
one as it uses the single characteristic time scale of the 
fluid evolution to define the velocities. Further the choice 
is essentially equivalent to that in cosmological simula- 
tions, as it means that the PS at small k (i.e. in the 
fluid limit) evolves exactly as in this context, growing 
in proportion to the square of the growing mode of lin- 
ear theory e'/ Tdyn 9 , and likewise the displacement and 
velocity fields in the fluid limit: 

u(R,i) = exp(t/r dyn )u(R, 0) (30) 
u(R,t) = exp(i/r dy „)u(R, 0) . (31) 



The more general form of the solution at leading order in this 
perturbative scheme n an in fact be derived using the formalism 
described here. See Till, 1131 for details. 

This is true because in the small k (fluid) limit only the curl-free 
component of the displacement field contributes to the density 
fluid fluctuations. This follows fr om t he relation Sp PS —V ■ u. 
which is valid in this limit. See JLl] , Il2j for the exact expressions 
and analysis of the fluid limit. 



III. PARTICLE COARSE-GRAININGS 

1. Defining a coarse-graining (CG) 

In this section we describe precisely how we construct 
the coarse-grained particle distributions we study. 

To coarse-grain a distribution of massive point parti- 
cles: 

• 1. We define a set of finite non-overlapping equal 
volume cells covering all space, i.e., a tiling of three 
dimensional space with equal volume tiles. We con- 
sider here, for simplicity, only the case that these 
cells are those of a simple cubic lattice, with lattice 
spacing which we denote £cg- Since the system 
studied here is evolved assuming periodicity in a 
cubic box of side L, we choose £cg — L/p where p 
is a positive integer. 

• 2. We ascribe one particle to each cell. The mass 
of this (point) particle is equal to the sum of the 
masses of the particles which are inside the given 
cell in the distribution being coarse-grained, and 
its position and velocity coincide with those of the 
center of mass and of these particles. 

We study here an initial particle distribution which is a 
SL or a SL evolved under gravity, in which all the points 
are ascribed equal mass. In general the point particles in 
the coarse-grained distribution just defined will not have 
the same mass since the number of points in the coarse- 
graining cells will not necessarily be equal. However in 
the cases which we consider in this paper, these parti- 
cle mass fluctuations introduced by the CG will be very 
small: (i) we will consider CG particle distributions at 
scales £cg which are large compared to the typical dis- 
placements of particles from their initial lattice position, 
and (ii) we will consider CG cells with 

l CG = a£,ae Z+ (32) 

which are perfectly aligned with the original lattice, and 
therefore each contain exactly the same number of points 
of the unperturbed lattice. To be more precise in quanti- 
fying the importance of these mass fluctuations, it will be 
useful to make the following distinction: one can consider 
the above described CG as a "Eulerian" CG, in contrast 
to a "Lagrangian" CG. In the latter case we modify step 
2 in the above definition to read: 

• 2'. We ascribe one particle to each cell. The mass, 
position and velocity, of this particle are those of 
the centre of mass of the points which were, before 
the application of the initial random displacements, 
on lattice sites inside the given cell. 

In this case — with CG cells aligned with the original 
lattice with £qg as in Eq. p2] ) — the CG particle distri- 
bution retains evidently the equal mass property of the 
original distribution. 
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2. Some properties of coarse- grainings 

A few basic results which follow trivially from these 
definitions are important to note: 

• A Lagrangian CG of a SL is itself an SL. 

• A Lagrangian CG of a SL with Gaussian PDF is 
itself a SL with Gaussian PDF. 

• An Eulerian CG of a SL converges, in the limit of 
small displacements, to the Lagrangian CG on the 
same grid. 

The first result is very simple to see. The position of 
the particles in the Lagrangian CG can be written as 

x/ = Rj + u(Rj) , (33) 

where 

1 (/) 1 (/) 

R ' = ^E R * > u ( R ') = ^E u ( R *) ( 34 ) 

i i 

and the sums are over the a 3 points i whose original 
lattice sites are inside the 7-th cell of the CG lattice. Rj 
is simply the center of the new CG lattice cell, and, since 
the u(R,-) are uncorrelated with variance A 2 , we have 

(u(Rj) • u(Rj)) - S K {I, J)^ = S K (I, J)A 2 CG . (35) 

a 

Thus the CG particle distribution is a SL, with variance 
reduced by a 3 . The PDF of the displacements in the SL 
can be calculated given that of the original particle dis- 
tribution, but in general has a different functional form. 
As noted above, however, for the case of a Gaussian the 
PDF is the unchanged, because of the fundamental prop- 
erty of stability of a Gaussian (see, e.g., Q). For the 
particle distributions we considered here as IC — SL with 
small amplitude Gaussian displacements — the Eulerian 
CG is practically equivalent to the Lagrangian one. 

We note one further important property of the CG 
particle distributions defined in this way: the PS of mass 
density fluctuations in the original particle distribution 
and that in the CG constructed from it converge at scales 
larger than ~ 2-7T /£cg- This follows from a well-known 
argument due to Zeldovich |^8[ [2^, J23|. This argument 10 
states that the perturbations to a mass distribution intro- 
duced by moving matter around on a finite scale, while 
preserving locally the center of mass, lead to a modifica- 
tion to the PS at small k (i.e. smaller than the inverse of 



We use the word "argument" because the heuristic derivation 
given by Zeldovich — and other authors in the cosmological lit- 
erature — fall well short of a proof. See [124] for a discussion 
of this issue and some explicit constructions which illustrate the 
result. 



Name 


N l/3 


L 


I 


A 


S 


m/ 17164 


SL64 


64 


1 


0.015625 


0.00195 


0.125 


1 


CG32 


32 


1 


0.03125 


0.00069 


0.0221 


8 


CG16o 


16 


1 


0.0625 


0.00024 


0.0039 


64 



TABLE I: Details of the three sets of initial conditions stud- 
ied in this section. SL64 is the original (or "mother" ) particle 
distribution and CG32o and CGI60 are the two "daughter" 
particle distributions obtained by CG it as described in the 
text. N is the number of particles and L is the box size (ar- 
bitrarily set to unity). The mass m of the particles is given 
in units of that in SL64, i.e., m^. The mass density is of 
course the same by construction in all the particle distribu- 
tions, po = Nm/L 3 = 1. 

the characteristic length scale) which is proportional to 
fc 4 . Since, as we have seen above, a SL has a PS which is 
proportional to k 2 at small k, the result follows. 

For the case of a CG on a SL, we can in fact derive this 
result directly from the expressions Eqs. ( |l2|) and Jl3| ) 
given above: for £cg = we have that no — > no/c? 
and, we have seen, A 2 — > A 2 jo? . Thus the leading k 2 
term of the PS, given in Eq. (|l3|), is invariant under the 
coarse- graining, while the next term, proportional to fc 4 , 
is not. 

This convergence of the PS of the particle distribution 
and its CG quantifies the sense in which the CG leaves 
invariant the large scale properties of the particle distri- 
bution. Since a continuous fluid-like description of such 
systems and their dynamics is expected to be valid for 
precisely such scales, one would anticipate that this de- 
scription should coincide for both systems. Indeed this 
is one of the underlying motivations for our study of this 
kind of CG. 



IV. COARSE-GRAINING OF THE INITIAL 
CONDITIONS 

In this section we study the evolution from a SL initial 
condition, and from CGs defined on this initial particle 
distribution. We first present and analyze our numeri- 
cal results, and then analyze what can be learned from 
them about the non- linear clustering in these systems. 
Specifically we will see that the non-linear two-point cor- 
relation function, which develops in a given such sys- 
tem, can be well reproduced, over a significant range of 
time and spatial scale, in a CG of the system evolving 
in a regime in which the relevant dynamics are driven by 
two body interactions of NN pairs. There is, however, 
a small discrepancy in the amplitude of the correlation 
function which we show can be quantitatively understood 
using PLT (and which is an effect of discreteness, i.e., due 
to the different particle density in the simulations). We 
then show that these, and other, differences between the 
mother (i.e. the original particle system) and daughter 
(i.e. the CG particle distribution) systems will increase, 
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without limit, as the parameter a = Icg/& does. This 
then leads us naturally to the CG presented in the fol- 
lowing section. 

As our original ("mother") system we take here, and in 
the next section, the following IC: a simple cubic lattice 
with N = 64 3 particles, to which we apply a shuffling 
with Gaussian PDF and normalized shuffling parame- 
ter 8 = 0.125, i.e., as given by Eqs. @ and (||) with 
A = 0.125£ (and I = L/N 1 / 3 ). We will call this initial 
condition SL64. 

We consider then also two other IC which are La- 
grangian CG as defined above: CG32o and CG16o, de- 
fined respectively on the 32 2 and 16 3 sub-lattice of the 
original lattice (i.e. with a — 2 and a = 4) respec- 
tively. Since 8 is small this CG is, as we have noted above, 
in practice equivalent to the Eulerian CG 11 . Given that 
A 2 /a 3 , we have 



a£ and & 2 CG 



^CG 

ice 



v5/2 



(36) 



which explains the values of <5 which appear in Table Q. 
Further, in order to keep the mass density po fixed, we 
take 



mcG = —? 

or 



(37) 



where m, the mass of particles in the original SL, is (ar- 
bitrarily) normalized to unity in SL64. 

The initial velocities of SL64 are fixed using Eq. (29) 



In the two CGs they are fixed, as described above, by the 
condition that the CG points are ascribed the velocity of 
the center of mass of the set of points they represent. It is 
evident that this is in fact equivalent (for the Lagrangian 
CG) to also using Eq. (|29| ) on the CG points directly. 



A. Results of numerical simulations 

The numerical code used for our simulations is, as in 
§ and Q, the publically available GADGET code ||. 
It uses a tree algorithm for the calculation of the gravi- 
tational force. The singularity at r = in the latter is 
regularized, for numerical efficiency, at very small scales, 
i.e., well below the initial inter-particle spacing £. The 
simulations here (as in Q ) use a "smoothing parameter" 
e = 1/15 ■ 1/64 fa 0.0010 in the units used in Table @, 
i.e., about one fifteenth of the smallest initial average 
inter-particle distance in our simulations 12 . We assume 



here, as in pi and indeed generically in all such cosmolog- 
ical simulations, that the results we consider as physical, 
starting from separations a few times larger than e, are 
essentially independent of this choice of e. As in Q we 
have tested this assumption, and found it to be good, by 
resimulating identical initial conditions changing only e 
to considerably smaller values. 

In Fig. ^ are shown snapshots of two-dimensional pro- 
jections of the simulations starting from SL64, CG32o 
and CG16o- Here, and in the rest of this paper, time 
is given in units of dynamical time T^yn- We stop our 
analysis at t = 8 as for longer times finite size effects be- 
gin to be important as the size of the structures formed 
approach the size of the box. As time evolves structures 
first form in each simulation at small scales — well below 
the initial grid size £ — and then subsequently at larger 
scales. Clustering develops first in SL64, and only later in 
its CG. While no similarity in the fluctuations is visually 
discernible when the particle distributions are grid-like, 
in the evolved configurations one sees clearly that the 
structures which formed in the CGs trace approximately 
those in the "mother" distribution. 

In Fig. ^| we show the evolution of the two-point corre- 
lation function £(r) at the times indicated. These plots 
reflect in a quantitative manner the information gleaned 
from the visual inspection of the snapshots in the pre- 
vious figure. The three systems show initially no dis- 
cernible similarity: the real space correlation function is 
in each case dominated by the features of the different 
underlying lattice. Strong clustering (with £ 3> 1) devel- 
ops first at small scales in the SL64 simulation, and then 
subsequently at larger scales. When the coarse-grained 
systems start to evolve a clustered regime, their corre- 
lation functions start rapidly to resemble strongly those 
in the "mother" system. We note, however, that there 
is a transient time in which the coarser system has a 
correlation function resembling that of the "mother sys- 
tem", but of a slightly lower amplitude. This "lag" of 
the coarser system is here manifest in the plot at t = 6. 
We will return to this point below. 

In Fig. U is shown the evolution of the PS in the three 
simulations at different times steps, with k in units of 
the Nyquist frequency — tt/£ of SL64. The three 
curves at t = correspond to the theoretical PS given in 
Eq. (pi), with the appropriate Gaussian form p(k) and 
value of the shuffling variance. The agreement of the PS 
at small k discussed above is clearly seen (and one can 
verify the accuracy of the k 2 behavior). The structure of 
peaks, which is different in each initial condition, comes 
from the second term in Eq. (O) 13 . Note that while the 



11 The probability p ou t that any given particle of the initial N = 
64 3 lattice is displaced by its shuffling so that it falls out- 
side the CG cell of the 32 3 CG can be written as p ou t = 
1 - [•\/ 3 / 2lr Jfi/M e~ 3x2 / 2 dx] 3 . For 6 = 0.125 this corresponds 
to p out ~ lO" 1 !. Thus Np ou t < 1. 

12 For this code |25| this means that the force is the exact grav- 
itational force for separations greater than 2e. For the precise 
functional form of the regularization below this scale, see |25|. 



13 Note that not all the Bragg peaks in the PS [cf. second term in 
Eq. ( |ll| )] are visible in the plot.) Indeed the first Bragg peak in 
SL64 is at k = 2kpf, while the first one in CG32o and CGI60 are 
at k = 4fcjv and k = 8fcjv respectively. The reason is that the PS 
is calculated using a sampling of fc, and the peaks are so narrow 
that they are missed. 
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FIG. 1: Snapshots at times (top line), 3, 6 and 8 (bottom line) of SL64 (left column), CG32o and CG16o (right column). 
What is shown is a projection on the x— y plane of a slice of thickness 0.6 along the z axis. 
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FIG. 2: Evolution of £(r) in SL64 and in CG32o and CGI60. The times are given in the plots. 



differences between the distributions at the initial time 
manifest themselves in the PS only at large k, we have 
seen that the initial correlation functions are very differ- 
ent at all scales 14 . At t = 3 this peak structure is still 
visible only in the CGI60 simulation, corresponding to 
the fact that the lattice structure is still present in this 
simulation which has evolved little non-linear clustering 
at this time. The relative difference between two other 
simulations is now less than at the initial time, show- 
ing again that once structures develop in the "daughter" 
particle distribution they trace those in the original dis- 
tribution quite accurately. At later times this behavior 
is also seen for the CGI60 simulation, as it "catches up" 



with the other ones. 

These observations are, to a first approximation, very 
much in line with the qualitative picture of the evolution 
of clustering widely accepted in the theory of structure 
formation perturbations in cosmology (see, e.g., |Q, or 
p7f for a concise summary): non-linear gravitational evo- 
lution transfers power from larger to smaller scales (by 
collapse) so efficiently that, for initial PS with n < 4 
(as is the case here, n = 2), the clustering amplitudes 
in the non-linear regime at a given time depend essen- 
tially only on the initial power in modes corresponding to 
larger scales. Thus what is relevant to recover the non- 
linear evolution of the system is to include this initial 
power at larger scales 15 . The coarse-grainings we have 



This "localisation" of the discrete characteristics of the distri- 
butions in reciprocal space, and "delocalisation" in real space, is 
also a feature of the initial conditions of cosmological simulations. 
Sec [Eql for a detailed discussion. 



See, e.g., |28[ for a, numerical study of this issue in cosmological 
simulation, and [M for a widely used phenomenological ansatz 
for "reconstruction" of the non-linear PS from the initial condi- 
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£/fc N (64 J ) 

FIG. 3: Evolution of the PS in SL64, CG32 and CG16 . 
Times are as indicated in the legend. The ir-axis is normal- 
ized to fcjy in SL64. At each time the different systems can 
be most easily distinguished by looking at the PS at large k, 
which approaches cowards 1/uq. The peaks which are evi- 
dent in in the PS are the Bragg peaks of the lattice, which 
progressively disappear as the lattice structure is washed out 
by the evolution. 



applied conserve, by construction, the initial power at 
larger scales, and thus can begin to reproduce clustering 
correctly at scales smaller than the CG scale once the 
smallest scales with the "correct" power go non-linear. 
More qualitatively, the approximate criterion used in cos- 
mology to determine when non-linear evolution should 
set in, adapted to the present case of a static universe, is 



k 3 P(k)e 



2t/Tdyn 



1 



(38) 



i.e., the fluctuations of wavenumber k satisfying this rela- 
tion at a given time are those "going non-linear" accord- 
ing to the linearized fluid theory at that time 16 . Assum- 
ing that the coarse-graining procedure leaves only initial 
fluctuations in reciprocal space below the Nyquist fre- 
quency, the largest k available to go linear, which is there- 
fore the first one to do so (since P(k) oc k n with n > —3), 
is the Nyquist frequency, = 2it/£cg- The latter de- 
creases, and therefore non-linearity sets in later, as Iqg 
increases. To determine the time scale this argument 
predicts, it is convenient to note that, using Eqs. Jl3| ) 
and (pf^), the criterion Eq. ( pq ) with k = kjy, may be 
rewritten 



Thus the time scale predicted for the first non-linear 
structures is just the time of fall of neighbouring particles 
on one another, calculated in the fluid limit (i.e. using the 



2t/r dyn 



1 . 



(39) 



ZA as in Eq. 27 ). We will return to the question of the 
accuracy of this prediction below, and explain that it will 
be a good approximation if discreteness effects up to this 
time (i.e. when particles first approach one another) are 
small. For the numerical simulations here this is indeed 
a reasonable first approximation, because the duration 
of this period is not too long (at most a few dynamical 
times). 

To the extent that the evolution of the measured two 
point quantities agrees (after some time), one can con- 
clude that this evolution can depend only on quantities 
which are the same in all the systems. Specifically, as 
we have just discussed, it can (and does) depend on 
the amplitude of the initial fluctuations at the relevant 
larger scales, but it does not depend on the particle den- 
sity. From this latter fact one might be tempted to con- 
clude that the dynamics of the system should be well 
described by one in which the density fields, for example, 
are smooth functions, i.e., by a description in which the 
discrete nature of the system is of no relevance 17 . Inde- 
pendence of measurable quantities of the particle density 
does not, however, establish such a conclusion. Indeed 
we will see now that, on the contrary, a description of 
the dynamics in a manifestly particle-like (i.e. discrete) 
way captures essential aspects of the evolution. 

We consider now more closely both the differences and 
similarities between the evolved simulations. We note in 
particular that the CG systems trace very well the two 
point correlations in the original system, already at a 
time when a very simple model for the development of 
these correlations, in terms of two body interactions, is 
valid. There is, however, a small discrepancy in the am- 
plitudes of the correlation functions, which we explain. 
This leads us naturally to consider the extrapolation of 
these results to the case that the CG scale becomes very 
large compared to the inter-particle distance in the orig- 
inal particle distribution. 



B. Non-linear correlations at early times 

In the results above we have seen clear evidence for the 
convergence of the correlation properties of the evolved 
coarse-grained IC toward those of the evolved "mother" 
particle distribution. There is, however, a range of time 
scales in each case in which there are already significant 
non-linear correlations in the CG simulation, but the 
agreement, notably in the amplitude of the correlation 



tions. 

The quantity k s P(k) is equal, up to a numerical factor, to the 
variance of mass fluctuations averaged over a spherical Gaussian 
window function of radius 2w/k cf . Jli| . We have adapted the 
criterion here to the case of a non-expanding universe, by using 
the linear theory growth factor appropriate for this case. 



We have in mind, specifically, the description of the dynamics by 
a Vlasov-Poisson system of equations. In cosmology the aim of 
N-body simulations is to reproduce the evolution of this system 
of equations. We will discuss this question in greater detail in 
the conclusions section. 
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FIG. 4: Reconstruction of using the NN distribution w(r) 
in Eq. (p"(i|). At time t — 2, it is applied to SL64: the thick line 
represents the directly measured £(r), while the thin line is the 
reconstruction with u>(r). The chosen times t = 4 and t = 6 
for the other two plots correspond to times a little before the 
correlation function of CG32o and CG16o catch up with SL64. 
In the two cases, we have plotted f(r) in SL64 (thick lines), 
£(r) in the CG system (dashed line) and its reconstruction 
with U)(r) (thin continuous line). This allows one to see that at 
these times, when the first non-linear correlations are created 
in all these systems, they are well accounted for purely by 
correlation of NN pairs, which develop through the interaction 
of such pairs. For clarity, the data at the two later times 
have been rescaled along the x-axis (by a factor 50 and 50 2 
respectively) . 



function, is not so good. This discrepancy can be un- 
derstood by analyzing the dynamics giving rise to these 
correlations at early times. 

We have shown in [|ll| that the evolution at early times 
— when the displacements are small compared to the lat- 
tice spacing — is very well described by th e per turbative 
approach (PLT) described briefly in Sect. 



[ID 



above. It 

breaks down when pairs of particles start to approach 
one another. In the immediately subsequent phase strong 
correlations (with £ ^> 1) begin to develop at small scales, 
and it is natural to hypothesize that an important role 
is played at this time by the interaction of particles with 
their NN. A simple way to test for this is, as explained in 
Sect. [IB above, is to see whether the relation Eq. jlO| ) 
holds. To do this we "reconstruct" the two-point corre- 
lation function £(r) from the measured NN distribution 
using this relation, and compare with the measured £(r) 
at each time. The results are shown in Fig. ^. Inspection 
shows that in each of the simulations the relation holds to 
a very good approximation. The times at which we have 
chosen to plot the comparison in each case is the largest 
time at which the agreement is observed to be good. This 
time increases as the IC becomes coarser because, as dis- 
cussed above, the phase in which PLT is valid increases 
as the CG scale does. Therefore the time of switching on 
of binary collapse becomes correspondingly larger. 



The CG systems thus develop significant non-linear 
correlation in the phase in which the relation Eq. jl0| ) 
holds. We infer that the correlation measured is es- 
sentially the result of two body correlations developing 
through the interaction of particles with their NN. What 
is very striking, as it has been pointed out in jl], [l5| , 
is that the two-point correlation function traces very 
well, at least in its shape, that in the mother simula- 
tion, even though in the latter the dynamics are at this 
time no longer NN dominated. At these times in fact 
the mother simulation is evolving in the self-similar man- 
ner described in Q, for which the temporal scaling [i.e. 
the dynamical exponent given in Eq. ([15|)] can be deter- 
mined from a fluid-like description of the system 18 . Thus 
the complex many body dynamics of the mother system 
(whether mean field or not) is well approximated, at the 
relevant scales, as particle-like "blobs" each moving un- 
der the effect of the single nearest "blob" . 



C. Extrapolation to large IcG/i 

There is, however, as we have noted also in the dis- 
cussion of Fig. ||, a small (but discernible) discrepancy 
between the amplitude of the correlation function in the 
CG system and that in the original one. The CG systems 
thus appear to lag behind for a time, before "catching 
up" subsequently with the mother system. Thus the dy- 
namics of the formation of the correlation in the mother 
system, in the asymptotic self-similar regime, appears to 
be very well represented by that of CG particles, but 
these particles fall toward one another by NN interac- 
tions slightly slower than the particle-like "blobs" of the 
original particle distribution. 

This difference can in fact be understood from the 



perturbative approach (PLT) of 1 1 1 , which we outlined 
above in Sec. II D. For most modes of the displacement 



field applied to a simple cubic lattice, the largest eigen- 
values characterizing growth are in fact smaller than the 
(single) exponent characterizing the fluid limit, given by 
the ZA as in Eq. ( |30[ ) for our initial conditions. Thus the 
evolution of displacements is, typically, slower in the dis- 
crete particle distribution (described accurately by PLT) 
than in the fluid limit 19 . Further the discrepancy in 
any given quantity in the two treatments becomes more 
important as the time increases (so long as PLT itself 
remains valid). Specifically, for example, the criterion 
Eq. (pt), derived above in the fluid (ZA) limit, for the 



18 We recall that the scaling is inferred from that of the large scale 
nctuations in the model treated in the linearized fluid limit, cf. 



19 We write "typically" because there are a very small fraction of 
modes, in a simple cubic lattice, which grow slightly faster than 
in the fluid limit. Given their small number, and the small differ- 
ence in the growth rate, they are significant in the average only 
at very long times. 
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time of onset of the phase of development of strong non- 
linear correlations (dominated, at first, as we have seen 
by NN interactions) can be calculated more accurately 
using PLT. The result can be written schematically as 



J CG C 



Jl-r,)2t/r dy „ 



1 . 



(40) 



where i] is different from zero, and typically positive 20 . 
Thus the discrepancy between the fluid prediction of 
Eq. (§|) and that of PLT Eq. @ will increase as 8 CG 
decreases, i.e., as a increases [cf.Eq. (36)]. 



It follows that if we extrapolate the coarse-graining to 
arbitrarily large a, an arbitrarily large difference in the 
evolution of the CG daughter and that of the original 
distribution will result 21 . Such behavior could, in prici- 
ple, be observed, of course, in numerical simulations like 
those performed here, if the initial 5 is sufficiently small. 
In practice this is not numerically feasible as it requires 
unattainable accuracy on the force 22 . We note also that 
the small time-lag effect which we have identified here in 
our simulations was not seen in the results of [p"5| , be- 
cause the SL simulations considered started from larger 
values of 5 (larger, specifically than that in CG16o, cf. 
Table 0). 



V. COARSE GRAINING AT A FINITE TIME 

We have seen in the previous section that the CG on a 
grid of the SL initial conditions, i.e., of the system at the 
initial time, leads to an evolution which diverges further 
and further from that of the original system as the ratio 
of the CG scale to the scale of the original grid increases. 
The reason is that the CG system of particles evolves for 
a longer time in the regime in which discreteness effects 
induced by the lattice accumulate. 

In this section we define and study a different CG 
which, by construction, avoids this problem: instead of 
coarse-graining the initial condition we coarse-grain the 



20 In [tlj the full gravity, PLT and fluid evolution is calculated 
for tne average relative displacement between particles initially 
on neighboring sites of a SL configuration, and the "lag" effect, 
described approximately by Eq. (ftCl) with r\ > 0, is manifest (see 
e.g. Fig. 8 of At very long times the sign of the effective 
rj will change as the few modes which grow more rapidly than in 
the fluid (see preceeding footnote) come to dominate. 

21 Indeed for asymptotically large a one will obtain completely dif- 
ferent qualitative behaviours of the two systems, as the small 
number of eigenmodes which grow faster than in the fluid limit 
will dominate. The eigenmode with the very largest eigenvalue 
(about ten percent larger than the fluid one) corresponds to ad- 
jacent pairs of planes of the lattice falling toward one another 
wM- In the limit that 5 is sufficiently small that such a mode 
can collapse before any other one, there will be an imprinted long 
range order in the structures formed which is sheerly an artifact 
of the CG lattice. 

22 In [|llj], however, related oscillating modes of the lattice have 
been simulated using a simplified code exploiting the symmetries 
of certain modes. 



original system when it has evolved for a finite time. The 
idea is simple: we choose the time at which we coarse- 
grain the system at a given scale by a prescription which 
means that, approximately (in the sense discussed in 
what follows), the parameter 5cg (i-e. the one point 
variance of the displacements of the CG particles off their 
lattice) is equal to that of the original particle distribu- 
tion. The CG system will therefore spend approximately 
the same amount of time in the regime of validity of PLT 
as the original system. 



A. Definition of coarse-graining at a finite time 

As in the previous section we consider coarse-graining 
a particle distribution with lattice spacing I using a CG 
grid which is also a lattice with lattice spacing 1 C g = a ^ 
where a is an integer. The prescription we use for the 
time tec at which we perform the coarse-graining is 



tec = Tdyn -In a. 



(41) 



The reason for the choice Eq. ( f4l| ) is simple : if the 
Zeldovich approximation is assumed to be exact at the 
scale of the coarse-graining, the CG system is then ex- 
actly a SL with the same value of the normalized shuffling 
parameter S. Thus the CG at time tec is related to the 
original system at t = by a simple rescaling of all length 
scales 23 . Let us see that this is the case. 

When the ZA is valid, in the form as given in Eq. j27|), 
it follows that all particle displacements are simply am- 
plified in time by the factor e*/ Tdyn . Thus a SL remains a 
SL, with a growing variance of its random uncorrelated 
displacements. If we perform a Lagrangian CG at time t 
we therefore obtain a SL with normalized shuffling Scg 
given as Eq. (|3tj), i.e., 



&cg = 



S(t) _ J(0)e*/ T ^- 



,5/2 



JS/2 



(42) 



Taking the prescription Eq. ( fi"l"| ) for the time at which 
the CG is performed, we thus have Scg = 8(0). Further 
we recall that in the ZA the velocities scale in proportion 
to the displacement [cf. Eqs. @ and Eq. (§|)]. Thus 
the velocities in the CG particle distribution also simply 
rescale those in the original distribution in proportion to 
a. The CG system is thus (statistically) identical to the 
original system up to a rescaling of spatial variables of 
the original system x — ► ax 24 . 

In reality of course the ZA is not exact. Above we have 
assumed in fact that it is exact for the evolution of each 
particle of the mother distribution until an arbitrary time 



23 This is for the infinite lattice and exact gravity. It is thus implicit 
here that we consider the regime where finite box size effects, and 
effects depending on the finite softening e may be neglected. 

24 We assume as the previous section that we have a Gaussian PDF. 
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(since tec increases arbitrarily as Icg does) . Clearly this 
is not even a good approximation in general. We need, 
however, only to make a much weaker assumption to re- 
cover the result (that the finite-time CG is a rescaling of 
the mother SL): all we require in fact is that the center of 
mass of the CG cell moves according to the ZA until tcG- 
Thus, specifically, we require only that the ZA be a good 
approximation when we neglect the gravitational forces 
acting between particles within the same CG cell (since 
such forces do not affect the center of mass motion). 

We have seen that the ZA can be derived as an appro- 
priate limit of the full particle evolution using PLT. This 
latter treatment is of course only valid until the time at 
which particles approach one another for the first time. 
We can, however, modify it to understand the nature of 
the approximation required here. To do so, we break the 
force on a particle into two parts: that due to particles 
in a region of size n <C ice about it, and the residual 
force from particles outside this region. The local part 
of the force will have a negligible effect on the motion of 
the center of mass of the regions of size Icg ^ V> and so 
we can neglect it. The second piece of the force can be 
linearized in the relative displacements as in PLT, which 
now remains until particles move a distance of order 77. 
Because the dynamical matrix has been modified only 
up to a finite scale, the diagonalisation which gives the 
eigenvectors and eigenvalues for the displacement field 
will converge to the same behavior at small k as in the 
full PLT. The ZA will therefore still be a valid approxima- 
tion in this limit, with corrections at any finite k which 
can in principle be calculated. Thus, in summary, we 
expect the ZA to be a good approximation as required, 
if the scale of non-linearity is smaller than £cg at time 
tec 25 - We also expect, as in the previous case of CG 
at the initial time, that, at a given length scale, the de- 
viations from ZA will diverge as the time for evolution 
with this "modified PLT" description does. We there- 
fore expect the result to break down if the initial S in 
the mother particle distribution (and therefore that in 
the CG distribution also, since they are equal) becomes 
itself arbitrarily small. 

The CG at a finite time, defined as we have just dis- 
cussed, is useful with respect to understanding the "self- 
similar" properties observed in the evolution of a SL, ex- 
tensively discussed in [|l] (and summarized in Sect. II C 
above). In the case that the CG is just a rescaling of 
the initial SL (i.e. in the approximation that the ZA is 
exact), it follows that 

Mr, t CG + 1) = £ (£, t) = C(re- 2tcG/5Td *'\t) (43) 

where £cg is the correlation function measured in the 



evolved CG system and £(r) that in the original one. 
Now, if the CG system indeed reproduces the evolution 
of the original system, i.e., ^cg(jlJ) = £,( r , t) for t > tec, 
the self-similarity relation Eq. (^5|) follows. 



B. Numerical simulations 

Our original ("mother") system is the same SL64 ini- 
tial condition considered in the previous section. Our 
coarse-grainings are now performed at the same length 
scales, i.e., a = 2 and a = 4 = 2 2 . Following the pre- 
scription just given we therefore coarse-grain the original 
particle distribution evolved to times tec and 2tcc re- 
spectively, where 



tcG 



-(ln2)T dyn 



1.73r d 



yn ■ 



(44) 



25 This will necessarily be the case here as the S of the original 
particle distribution, which controls the time at which any scale 
goes non-linear, has been chosen small compared to unity (but 
not too small, as discussed above). 



Unlike the case of the CG at t = 0, the Eulerian and 
Lagrangian CGs are in each case not the same. We thus 
define for each case two CGs so that we have four in to- 
tal CG32 t , CG32J, CG16 t and CG16' t . The unprimed 
acronyms refer to the Eulerian (non-equal mass) simula- 
tions, the primed ones to the Lagrangian (equal mass) 
simulations. Since the IC are actually such that particles 
have moved only a small distance compared to Icg at 
the time of the corresponding coarse-graining, the differ- 
ences should not be very important, as they come only 
from particles close to the surface of the CG cell being 
re-attributed to neighboring cells. 

In Fig. ^| we give a schematic representation of this 
CG scheme. In Fig. || are shown snapshots in time of 
the evolved SL64, CG32 and CG32J systems starting 
from a time t — 2 just a little after tec ~ 1-73 at which 
the latter distribution is defined. CG32o and CG32J are 
clearly very similar, but visually one can discern that 
CG32J traces slightly better than CG32 the morphology 
and dumpiness of the structures in the mother particle 
distribution. 

In Fig. [7] is shown the evolution of the two-point cor- 
relation function in SL64, CG32J and CG16J. We do 
not show here results for the Eulerian CGs because the 
(number-number) correlation function we can estimate 
straightforwardly is not equal in this case to the density- 
density correlation function (which is the quantity we 
are interested in). As in the previous section we see that 
the CG particle distributions begin to develop non-linear 
correlations later than in the mother distribution, but 
that the relative difference between the correlations in 
the two particle distributions becomes negligible subse- 
quently. In Fig. H is shown the temporal evolution of 
the PS in the different particle distributions. In this case 
we show both the Lagrangian and Eulerian CG distri- 
butions as our estimator of the PS is indeed that of the 
mass density fluctuations [i.e. it takes into account the 
different masses of the particles, cf. Eq(^)]. We see again 
qualitatively very similar results to those in the previous 
section. The CGs, when they are initially defined, have 
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FIG. 5: Schematic representation of the "CG in time" exper- 
iment. The 64 3 , 32 3 and 16 3 systems are represented from 
top to bottom, and the time evolution is represented in steps 
of icG from left to right. The arrows marked "CG" indicate 
that a CG is done (note that in the experiment done, the 
CG for the 16 3 system is obtained from the 64 3 system and 
not from the 32 3 ). The double arrows marked "1. s." (large 
scales) indicate that the system should share the same fluc- 
tuations at sufficiently large scales. The dashed arrows and 
the dashed squares illustrate the "self-similarity" argument: 
a square defined (partially or wholly) by the dashed lines are 
supposed to be identical up to a rescaling of length scales to 
the system to which the dashed arrow points. 



small correlations which arc indistinguishable in the fig- 
ures from those in the mother particle distribution, but 
at larger k (of order the inverse of the CG scale Icg) they 
differ. We see also the small differences at these scales 
between the Eulerian and Lagrangian CG. As in the pre- 
vious section we see that the relative differences between 
the evolved particle distribution decreases in time, with a 
greater convergence as non-linearity develops at smaller 
scales. Each PS must converge to the corresponding Pois- 
son noise level at large fc, which is simply the inverse of 
the particle density. Outside this range the PS of each 
CG interpolates rapidly, in the evolved particle distribu- 
tions, to the PS of the mother distribution. 

As discussed above, the CG at a finite time we 
have defined is useful for understanding the self-similar 
properties of the asymptotic evolution from SL initial 
conditions Q . Indeed it has the property that, if the ZA 
is exact, the CG particle distribution defined in this way 
is identical (statistically) to the initial condition up to a 
rescaling of length scales. The "self-similarity" is then a 
direct consequence of the fact that the evolution before 
the onset of non-linearity on a given scale depends to a 
very good approximation only on the initial fluctuations 
at these scales, which themselves have (by construction) 
no intrinsic length scale, and the purely fluid evolution 
of these fluctuations described by the ZA. 

In this respect it is interesting to check directly to what 
extent the CG mass distribution actually is a rescaling 
of the mother distribution, i.e. to what extent the ZA 
is an accurate approximation to the motion of the cen- 



ter of masses at the scale £cg, by direct comparison of 
the particle distributions. Further it is interesting to see 
how well the relation Eq. ( p3| ) is obeyed. In Figs. |9| and 
|j"o| we show these comparisons. In the first the dimen- 
sionless combination no-P(fc) is plotted as a function of 
the k wavenumber, rescaled as required to the Nyquist 
frequency in each particle distribution (this corresponds 
to a rescaling of length scales by a). In Fig. [Hj like- 
wise the x axis has been rescaled in the same manner 26 . 
We see that the particle distributions are indeed initially 
very similar — except for some visible differences close 
to the Nyquist frequency, due to deviations from the ZA 
at these small scales — and in time become more similar 
as the evolution develops. 



VI. CONCLUSIONS AND DISCUSSION 

We have studied the evolution of an infinite SL of 
points under its self-gravity by comparing it with that 
of its particle CGs, using two different prescriptions for 
the latter. The primary aim is to get insight into the dy- 
namics of the non-linear gravitational clustering in this 
system, which has been shown in to be very similar 
qualitatively to that in simulations of structure formation 
in cosmology. We first summarize the main conclusions 
we draw from our study. We then discuss one specific 
question, which is important with respect to the theoreti- 
cal interpretation of these results: the definition of a limit 
in which a continuum theoretical description should be 
appropriate to fully understand the evolution (in the lin- 
ear and non-linear regime) of this system. This question 
is of particular relevance in the context of qualitatively 
similar simulations in cosmology, as such simulations are 
useful insofar as they reproduce the evolution of such a 
limit. How they do so is, however, currently poorly un- 
derstood. Finally we briefly mention some directions for 
further work. 



A. Summary of results 

Our principal results and conclusions are the following: 

• In the numerical simulations the coarse-grained 
particle distributions are observed to evolve to give, 
after a sufficient time, two-point correlation prop- 
erties which agree well, over the range of scales 
simulated, with those in the original distribution. 
Indeed both the original system and its coarse- 
grainings converge toward a simple dynamical scal- 
ing ("self-similar") behavior with the same ampli- 
tude. The characteristic time required for the CG 



See the footnote in discussion of Fig. ^ above for the explanation 
of why not all the predicted Bragg peaks are visible in the plot. 
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FIG. 6: Snapshots of the evolution of SL64 (left column), CG32o and CG32j (right column) at times 2,4,6 and 8. The points 
shown are obtained by projection on the x-y plane of a slice of thickness 0.4 along the z axis. 




FIG. 7: Evolution of £(r) in SL64 and the Lagrangian (equal mass) coarse-grainings CG32J and CG16(. The different plots 
correspond to the times 1, 2, 3 and 4 in units of tcG- 



system to begin to reproduce the clustering in the 
original particle distribution at scales below the CG 
scale increases as the latter scale does. As discussed 



in Sec. [VA these observations are all very much 
in line qualitatively, and to a first approximation 
(see below) quantitatively, with the qualitative pic- 
ture of the evolution of clustering widely accepted 
in cosmology: the CG distributions share the same 
fluctuations at large scales and it is these initial 
fluctuations alone, to a very good approximation, 
which determine the correlations which develop at 
smaller scales at later times. 



• Once particles begin to fall on one another in any 
of the SL simulations, there is a phase in which 
very significant non-linear correlations develop due 
to interactions between NN pairs of particles. Fur- 
ther we have seen, as in jjl that the form of 
the two-point correlation function which develops 



in this phase is very similar to that observed, in the 
same range of amplitude, in the asymptotic scaling 
regime at later times. Thus it appears that it is al- 
ways possible to choose a CG of the original system, 
which reproduces quite well the non-linear correla- 
tions in the original system with this "early time" , 
explicitly discrete, dynamics of "macro-particles" 
of the CG distribution. 

This provides a simple physical picture/dynamical 
model for the generation of the non-linear correla- 
tion function in the relevant range (up to an am- 
plitude of £ in the range of 10 — 10 2 ). We will 
describe more quantitatively in a forthcoming pa- 
per 1 30 the exact form of the non- linear correlation 
function given by this model. 

This finding is prima facie very different, in spirit 
and substance, to any existing explanations of the 
dynamics giving rise to non-linear correlations in 
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FIG. 8: Evolution of the PS in SL64 and the different CG in time (see the legend). The different plots correspond to the times 
1, 2, 3 and 4 in units of icG- 



N body simulations in cosmology. In this context 
theoretical modeling invariably assumes that the 
non-linear correlations observed in simulations in 
this range should be understood in the framework 
of a continuum Vlasov limit, in which a mean-field 
approximation of the gravitational field is appropri- 
ate. Indeed the fact that self-similarity is observed, 
with a behavior independent of the particle density, 
is usually taken as an indication that such a contin- 
uum description is appropriate. Our model is man- 
ifestly not of this type (while also consistent with 
the amplitudes of the correlation function being in- 
dependent of particle density). There is, of course, 
not necessarily a contradiction: it may be that our 
simple discrete model is simply a sufficiently good 
approximation to the continuum model 27 . We will 



return to a discussion of this point below. 

• A closer examination of the results shows that there 
are small, but observable, differences in amplitudes 
between the original SL and it CGs. These may 
be explained using a more accurate description of 
the early time evolution of the system, given by 
the perturbative PLT treatment of O, valid 
up to the onset of the NN interaction phase. In 
this framework the fluid limit criterion Eq. (|3^) for 
the end of this phase may be replaced by a more 
accurate criterion, taking into account the depen- 
dence of the rate of evolution of modes of the dis- 
placement field at small scales on particle density 
(and their wavelength). Generically the effect of 
discreteness (i.e. the difference with respect to the 
fluid limit) is to slow down the growth of modes 



In [Rlfl we have presented an analysis of the force distribution in 
cosmological simulations at the time of formation of the first non- 
linear correlations, showing that NN interactions also dominate 
in this case. Thus we do not expect that either the different 



initial conditions, or the fact that there is expansion of the spatial 
background, should qualitatively change the early time dynamics 
with respect to that of the system discussed here. 
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FIG. 9: Comparison of the PS in SL64 at times 0, 1, 2, 3 in units of tec with the PS in CG32J at times 1, 2, 3, 4 and the PS 
in CG16t at times 2, 3, 4, 5. We have plotted noP(k) as a function of k/kN, where fcjv is the Nyquist frequency for the given 
distribution. In the approximation that the ZA is valid at the scale of coarse-graining, all the plots should overlap. 



at smaller scales (comparable to the inverse of the 
inter-particle spacing). Thus the time of exit from 
the PLT phase is longer than given by Eq. ( p9| ) . In 
our numerical simulations a resultant "lagging" of 
the non-linear correlation function in the CG sim- 
ulation is just discernible at early times. These ef- 
fects are small, and visible only at very early times, 
because the value of S used, albeit small compared 
to unity, is not so small. In the limit that the 
normalized shuffling 5 becomes arbitrarily small, 
however, these differences between the discrete and 
fluid system can become arbitrarily large. 

• Since the coarse-grainings defined on the initial 
conditions give 5 oc (£/(cg) 5 ^ 2 , differences between 
the full PLT evolution and the fluid evolution will, 
for the reasons explained in the previous paragraph, 
diverge in the limit {I/Icg) — * 0. The coarse- 
grainings defined at finite times avoids this feature 
by construction: it ensures that the CG system has 
approximately the same S, so that the time spent 
in the phase in which discreteness effects accumu- 



late as described by PLT is identical to that in the 
mother particle distribution. Thus the use of a CG 
distribution, evolving in its "early" phase, to re- 
produce the non-linear correlations in the original 
distribution would be expected to work in this case 
even for arbitrarily large times (which require arbi- 
trarily large values of £cg)- 

• The origin of self-similarity observed in the evolu- 
tion of the system, in which the temporal evolution 
of the correlation function is equivalent to a rescal- 
ing of spatial variables, is nicely illustrated using 
the CG at finite times. Indeed insofar as the sys- 
tem, for its further evolution, can be approximated 
by its CG defined in this way, one obtains self- 
similarity. We have seen that there are thus two 
ingredients giving this behavior. Firstly the non- 
linear correlations which develop at a given scale 
depend on the initial conditions only through the 
fluctuations present initially at larger scales (i.e. 
the same as the first point above). Secondly these 
larger scales are well described at early times by the 
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FIG. 10: Comparison of £(r) in SL64 at times 0, 1,2,3 in units of tcG with £(r) in CG32J at times 1,2,3,4 and in CG16£ at 
times 2, 3, 4, 5. As in the previous figure, the approriate rescaling of lengths on the x-axis has been done so that these functions 
should overlap in the approximation that ZA is exact at (and above) the CG scale. 



Zeldovich approximation, i.e., by the purely fluid 
limit in which the evolution is independent of the 
particle density. 



B. The continuum limit 

As noted in the introduction, one of the motivations 
for studying coarse-grainings like those we have consid- 
ered here is that such a study may help to clarify what 
properties of the evolved SL may be understood as valid 
in a continuum limit: the idea is that such properties 
should be invariant under coarse-graining, i.e., invariant 
with respect to changes in the small scale discrete struc- 
ture of the systems. The fact that we have found that 
the evolution observed in the original system is not com- 
pletely invariant in this sense, and indeed that the dif- 
ferent coarse-graining procedures give different results, 
shows that this evolution is not fully representative of 
such a continuum limit. The question of such a limit is 
of relevance if one wants to determine whether, in par- 



ticular, the form of the correlation function in the self- 
similar should be derivable with such a continuum the- 
ory. As mentioned above, this is particularly important 
to understand in the context of simulations of this type 
in cosmology, as the goal of the N body method is to 
reproduce such a limit. 

To address this question we need to make more precise 
what we mean by "continuum limit" . There is in fact, of 
course, no unique way of defining such a limit. Implic- 
itly above we have been supposing a limit in which the 
inter-particle spacing t goes to zero in some appropriate 
manner, so that a hydrodynamic type description of the 
system becomes valid. More specifically, in cosmology 
and more generally for systems with long range interac- 
tions, one is interested in identifying the Vlasov limit, 
in which the system is described by a one particle phase 
space density obeying a Liouville equation in which the 
force is calculated in a self-consistent mean field approx- 
imation (see, e.g. (ll| 32 ). 

The SL system we have studied is characterized [[t| by 
a single dimensionless parameter: the normalized shuf- 
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fling parameter 5. Indeed the evolution we have studied 
is, as discussed in 0], independent of both the system 
size and force softening. Two SL with the same S, but 
different £, are then equivalent because of the scale free 
nature of gravity. To give meaning to a limit /->0we 
must, therefore, evidently introduce at least one addi- 
tional length scal e. W e note that in the discussion of the 
ZA above, in Sec. HE, we have derived it as a continuum 



limit by fixing the wave length of the mode considered 
and taking t — > 0. In this way we established that long 
wavelength modes of the displacement field behave as de- 
rived from a set of fluid equations for the self-gravitating 
fluid. The SL system as a whole, however, cannot have 
any continuum limit in the desired sense unless we intro- 
duce an additional length scale in the system. 

In fl33| it is shown formally that such a Vlasov descrip- 
tion is recovered in the limit that N — ► oo for a class 
of long-range interacting system of N particles, provided 
their interaction potentials are regulated at small sepa- 
rations. This suggests that one way in which we should 
recover such a continuum description is by introducing 
such a regularization, characterized by a scale e. Then 
one can take £ — > at fixed e (and then, in principle 
e -> 0). 

Following what was said above, we expect that in this 
limit it should be possible to define coarse-graining pro- 
cedures like we have done which leave the evolution of 
the system invariant. It is not difficult to see that this 
should indeed be the case. Let us assume that the force 
falls rapidly to zero at separations less than e. The PLT 
treatment can be applied and will give a spectrum of 
modes which are modified with respect to full gravity for 
fee > 1. In the limit £/e — > these modifications will 
become independent of I. Thus the evolution of the dis- 
placement field will be independent of I. As the CG de- 
fined are just a modification of this latter scale, they will 
have no effect on the evolution (provided £cg/z -C 1)- In 
practice this means that we can coarse-grain the system 
of finite £ (and t <§; e) up to any £cg, provided £qg "C e. 

An alternative way of defining the limit £ — > is by 
introducing a characteristic scale not in the force, but 
in the fluctuations in the system. The SL has no such 
scale, but it is easy to introduce one by defining instead 
the initial displacement field so that it has such a cut- 
off: we derive it as the sum of modes in reciprocal space, 
with all modes k > k c set equal to zero. We can then 
define a continuum limit for the system as k c £ — ► 0. It is 
not difficult to see that a similar reasoning to that used 



for the previous case leads to the conclusion that coarse- 
graining like those we have defined now change negligibly 
the evolution of the system provided k c lcG <C 1- 

We note that there is one intrinsic problem with the 
latter prescription: it works because the scale k c is the 
one characterizing the shortest wavelength fluctuation in 
the system. This is only true, however, at the initial 
time, as gravity will progressively (and very efficiently) 
create fluctuations at smaller scales. Thus the continuum 
limit will only be valid for some finite time. For the 
former case this will not be the case as the finite value 
of e will place a time-independent lower bound on the 
scales at which fluctuations develop. There is thus an 
intrinsic time independent separation between the scales 
of inhomogeneity and the particle scale. 

This discussion tells us then how we should extrap- 
olate our simulations toward a well defined continuum 
(Vlasov) theory to determine which results are valid in 
this limit: we should increase particle density until l<e. 
Manifestly this will "kill" the discrete particle dynamics 
which we have used in our model to explain the form 
of the non-linear two-point correlation function. If this 
form is indeed that of the continuum limit it should be 
observed in simulations with l<£ when non-linear cor- 
relations develop for r S> s. To perform simulations to 
test whether this is so requires that there be a reasonable 
range of scale r ^S> e ^S> £ well inside the simulation box, 
which implies a particle number much greater than that 
in the largest 64 3 simulation discussed here (Taking, e.g., 
I = 0.2s, e is already almost one tenth of the box size). 
With the particle numbers, in the range 10 s — 10 10 , now 
being accessed by several gravitational N body simula- 
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tion groups (see, e.g. 
are, however, in principle feasible and would yield clear 
conclusions on this very basic question. 
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